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on the free energy. It is illustrated by application to the low-concentration/high- 
concentration transition in globular proteins where a novel two-step mechanism is 
identified where the first step involves the formation of long-wavelength density fluc- 
tuations and the second step is the actual nucleation event occurring within the 
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I. INTRODUCTION 



Nucleation is a fundamental physical process of importance in fields as diverse as chem- 
istry, materials science, biology and cosmology. Our basic understanding of it goes back 
to Gibbs' discussion of the nucleation of a liquid droplet in a metastable vapor, see e.g. 
discussion in Ref. Q. The physics is governed by the fact that the excess free energy of a 
cluster of new phase relative to the mother phase has a negative contribution that scales 
as the volume and a positive contribution that scales as the surface area. The former is 
due to the fact that the bulk new phase has lower free energy than the bulk mother phase 
while the latter arises from the free energy cost of the interfacial region. When the cluster is 
small, the surface term dominates and the excess free energy of the cluster is positive mak- 
ing the cluster thermodynamically disfavored. When the cluster is large, the volume term 
dominates, the cluster has lower free energy than the vapor and so cluster growth is favored. 
Separating the two regimes is the metastable critical cluster. In Classical Nucleation Theory 
(CNT) it is assumed that the cluster is spherical, that its interior is in the bulk new-phase 
state and that the surface tension is the same as for a planar interface of the coexisting new 
and mother phases so that the free energy of the cluster as a function of its radius can be 
calculated giving a quantitative picture of homogeneous nucleation^. 

The dynamics of nucleation - the nucleation rate, the determination of the nucleation 
pathway and the nature of the post-critical growth of clusters - is a more involved sub- 
ject. In CNT, the starting point for understanding nucleation rates is the monomer attach- 
ment/detachment picture that leads to the Becker-Doring equations^. These are simply 
rate equations for number of clusters of size N based on the processes of monomer attach- 
ment and detachment. Following Frenkel, Zeldovich, and others^, these equations can be 
converted into a Fokker-Planck equation for the evolution of the distribution of clusters. 
All of these formalisms can be used to estimate nucleation rates. The physics enters via 
the monomer attachment and monomer detachment rates which in general depend on the 
thermodynamic driving force (the derivative of the CNT free energy with respect to cluster 
size) and the dynamics of the surrounding mother phase. In many cases, adequate approxi- 
mations are available that express the nucleation rate in terms of e _/3AQc where (3 = 1/ksT, 
ks is Boltzmann's constant, T is the temperature and AQ C is the excess free energy of the 
critical cluster. 
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Although physically appealing, this formalism suffers from a number of deficiencies. First 
and foremost is the well-known fact that the nucleation rates estimated from CNT are often 
many of orders of magnitude too large or too small^ (see also Ref. [l, Fig. 13.8). This is 
largely attributed to a poor estimate of the free energy barrier: more sophisticated methods 
allow for better estimates of the barriers. A common way to get improved estimates of 
the barriers is via classical Density Functional Theory (DFT)^^ (for early examples of this 
approach see Refs. SandS). However, the use of DFT raises new conceptual questions. The 
monomer attachment/detachment picture requires estimates of the rates for all cluster sizes 
which in turn means estimating "free energies" for all clusters. But, DFT is an equilibrium 
theory that can only describe - at best - stationary states on the free energy surface. This 
fact follows trivially since the procedure to determine free energies via DFT begins with 
the extremization of the free energy to get the structure of the cluster and naturally must 
result in a minimum, maximum or saddle point but it is also connected to the theorems that 
underlay DFT, which only attach significance to minima of the free energy functional*^. 
This therefore raises the question as to whether DFT can be used to determine the nucleation 
pathway - i.e. the sequence of necessarily nonequilibrium states leading from the pure 
unstable phase to the critical cluster. 

Putting aside the conceptual questions, there are currently two primary approaches to 
determining nucleation pathways using DFT. The first, and oldest, is to extremize the free 
energy under some sort of constraint that stabilizes non-critical clusters: in effect, imposing 
a stabilizing external fielcr^ii. The problem with this approach is that there are numerous 
possible ways to implement it. To take a simple example, one could minimize the free 
energy at fixed radius so that the pathway would be mapped out by varying the radius 
which acts as the reaction coordinate (or "order parameter"). However, one could equally 
well minimize at fixed number of molecules in the cluster, N, so that this becomes the 
order parameter. These do not necessarily give the same result since, e.g., the number of 
molecules can increase with the radius actually decreasing provided the interfacial width or 
interior density of the cluster increases enough to compensate; indeed, this actually occurs in 
some calculations^. Notwithstanding this problem, the constraint approach continues to be 
widely usecr^r— . A more recent alternative is to search for steepest descent paths on the free 
energy surface^ 1 ^ - — . However, here there is also an ambiguity as steepest descent requires 
a measure of distance in parameter space and there is a priori no unique prescription for 



this measure^. 

It is reasonable to suppose that all of the ambiguity underlying the DFT approaches 
is due to the fact that DFT is an equilibrium theory whereas nucleation is fundamentally 
a nonequilibrium, fluctuation- driven process. Viewed in this way, one can question other 
aspects of CNT as well, such as whether it is really necessary that the system even passes 
through the critical cluster since the free energy may not be the only factor controlling the 
nonequilibrium evolution. Similarly, one can also ask why the pre-critical process does not 
seem to depend on fluid-kinetics while these are crucial in driving post-critical growth. (This 
question has also been addressed recently by Peters^ using CNT free energy models). All 
of this suggests that one would like a unified description with the following properties: it is 
inherently nonequilibrium in nature, it is dynamical; it allows for the determination of the 
nucleation pathway from the dynamics, it allows for the determination of nucleation rates 
and it gives a description of post-critical growth. Of course, one would expect that it involves 
free energies in a natural way without assuming equilibrium and that under some limit the 
CNT would result thus also helping to clarify its applicability. Additionally, one would hope 
that such a formalism would be sufficiently flexible to be applied to more complex problems 
involving multiple order parameters and to determine the relative probability of completely 
different nucleation pathways (e.g. whether crystals form directly from solution or whether 
they form in two steps, passing first through a dense amorphous phased). 

The aim of this paper is to illustrate an approach that can satisfy these requirements. 
Attention will primarily be focused on nucleation from a vapor or liquid phase, but the basic 
ideas could be extended to other systems in obvious ways. The development is based on 
Brownian Dynamics where molecules move according to Newton's laws while being subject 
to a frictional force as well as fluctuating forces. This is a simple model for colloids and the 
important case of macromolecules in solution in which cases the friction and the fluctuations 
come from the bath/solvent and has the virtue of allowing for useful simplifications. As 
described in Section II, the fundamental level of description adopted for the description of 
nucleation is fluctuating hydrodynamics, in which the free energy enters the dynamics via 
a local equilibrium expression for the pressure, resulting in a set of stochastic differential 
equations (SDE's). The Most Likely Path (MLP) is introduced as a means of characterizing 
the nucleation pathway via a deterministic equation derived from the stochastic model. In 
Section III a reduced description in terms of order parameters is developed. It is shown 
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that in the weak-noise limit the MLP is simply gradient descent on the free energy surface 
with a metric uniquely determined by the dynamics and that it necessarily passes through 
the critical cluster thus making contact with the CNT picture. Finally, in Section IV three 
applications are discussed: nucleation of high concentration droplets from a low concentration 
solution of globular proteins; the reverse transition of high concentration solution to low 
concentration bubbles in an under-saturated solution and finally the specialization to a single 
order parameter so as to recover well-known results from CNT. For the first application 
- droplet nucleation - calculations lead to an unexpected but physically appealing result 
wherein it is found that the coupling of fluctuations, fluid dynamics and thermodynamics 
implies that nucleation is most likely to occur within regions of enhanced density due to 
long-wavelength density fluctuations. The final Section summarizes the results and discusses 
possibilities for further developments. A brief summary of some of these results has appeared 
previously^. 



II. THE MODEL 
A. Dynamics 

The system consists of a collection of particles with positions and momenta qj,pj inter- 
acting via a potential V. Additionally, the particles interact with a bath/solvent of small 
particles and this is described via a frictional drag and a fluctuating force, both characterized 
by the constant 7, 

<U = Pi/m, Pi = -7^- - 7Pi + e^2^mk B Tii (t) (1) 
where all the components of the force are independent, 

(S(t)£,-(0> = iM(*-0 (2) 

The coefficient of the noise term is determined by demanding that the system admits of 
an equilibrium Boltzmann distribution thus forcing a fluctuation-dissipation relation. Here 
and in the following, an additional parameter, e, is included to allow for the discussion of a 
weak- noise limit. The physicality of such a limit will be discussed below. 
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Defining the local density and momentum density respectively as 



r - q< 



(3) 




one sees that these satisfy the exact equations 
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with 



(F fat) F fa;t')} = 5(t-t')5(r- r') 1 



(5) 



Note that the actual form of the noise term that occurs in the momentum equation is 



However, it is well-known that a noise term of the form JB^^-(t) (with summation over 
repeated indices) can be replaced by another of the form without changing the 

statistical properties of the dependent variables provided that coefficients are related by 
BijBji = CuCji, i.e. the autocorrelation of the noise is unchanged^ 1 ^. 

There are two sources of stochasticity in this problem. The obvious source is the fluctuat- 
ing force that represents the bath or solvent. The second arises because one is not interested 
in following a single trajectory of the system but, as usual, want to consider an average over 
an ensemble of systems. In equilibrium, the ensemble is drawn from a known distribution 
(the canonical distribution). The distribution for systems out of equilibrium is generally 
time- dependent and one must specify the ensemble for a particular point in time, e.g. a dis- 
tribution of initial conditions. In the present case, because of the two sources of noise, one 
must carry out this procedure in two steps: first, the set of values of the fluctuating forces 
acting on the colloids at each instant in time is specified. Then, this fixed set of values for 
the noise simply acts as an external field coupled to the colloid dynamics and nonequilibrium 
statistical mechanics can be developed in the usual way as an average over an ensemble of 
initial conditions of the particles. Then, having developed whatever level of description is 




(6) 
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of interest, a second average is performed over the noise. As a result, one can go from the 
exact balance equations, Eq.(JlJ, to a description of ensemble-averaged quantities. 

To fix the notation, let T(t) = {<li,Pi}iLi represent the collection of positions and mo- 
menta of all iV large particles at time t and to be more explicit one can write r(t;r(0)) 
since the phase variables at time t are completely determined by their values at some earlier 
time. Similarly, one can write for the density, for example, p(r;t) = p (r; T(t; T(0))). The 
probability for the phase variables to have particular initial values is denoted / (r(0)) and 
the ensemble-averaged density is then 

p(r, t) = Jp (r; T(t; T(0))) / (r(0)) dT(0) (7) 
= J p(r;r)f(T;t)dT 

where in the second line the time-dependence has been moved to the distribution function^ 5 .. 
Notice that the ensemble-averaged density and momentum density are all conditional on 
the fluctuating forces which for now just act as external forces. At this point, the usual 
methods of nonequilibrium statistical mechanics can be used to develop a hydrodynamic 
description - either a mesoscopic fluctuating hydrodynamics or a macroscopic, deterministic 
hydrodynamics, the former having the form 



dt 



^-iV.j(r) = (8) 



at , V ■ j (r) j (r) lp (r) + pVp(r, t) + V • II (r) = - 7 j (r) + e^2m^k B Tp (r)£ (r; t) 

where p (r) and j (r) are the ensemble-averaged local density and momentum density and 
II(r, t) is the dissipative stress tensor which for the case of fluctuating hydrodynamics will 
have both a deterministic part and a fluctuating part. The other new term that appears 
here is p(r, t) which is the local pressure. Note that in the present case, the bath acts as a 
thermostat and so temperature fluctuations and heat transport are neglected, although in 
principle the model could be extended to include them. 

The final approximation involves the pressure and can be motivated in different ways but 
always amounts to an assumption of local equilibrium. At its most heuristic, one starts with 
the Gibbs-Duhem equation for a single-component system, Ndp = —SdT + Vdp. For an 
isothermal system, this gives dp = pdp where p is the density. Assuming local-equilibrium 
this is generalized to 

dp(r, t) = p(r, t)dp{r, t) = p(r, t)d^j^ (9) 



which is a widely used approximation (see, e.g., Ref. [l| Section 8.4). Substituting into the 
equations for the density and momentum gives 

^ + ^V-j(r) = (10) 
^ + V ■ j (r) j (r) /p (r) + p (r) Vfg + V • II (r) = - 7 j (r) + e^lm^Tp (fj£ (r; t) 

An alternative way to motivate this is to evaluate the pressure under the assumption of a 
local equilibrium (i.e. Maxwellian) velocity distribution and using the local-equilibrium ver- 
sion of the first YBG equation^ to rearrange the resulting virial expression for the pressure. 
The same system of equations has been developed starting with the Boltzmann equation^!. 
Chavanis^ has discussed the derivation of these equations for coarse-grained quantities (i.e. 
averaged over small volumes and time intervals) as well as their derivation based on the gen- 
eral Landau-Lifshitz fluctuation formalism^. Chavanis also discusses the difference between 
this model and the Cahn-Hillard dynamics, the main point being that the latter does not 
include the density in the noise amplitude and the free-energy term. However it is found, the 
most important assumption is clearly the introduction of the local equilibrium expression 
for the pressure since this is what brings the free energy into the picture. This assumption 
of local equilibrium is, on the one hand, heuristic and represents an uncontrolled approxi- 
mation. On the other hand, it is widely used in nonequilibrium statistical mechanics (e.g. a 
similar assumption underlies the generalized Enskog equation^) and similar approximations 
are even used for dissipative granular systems. Thus, while it cannot be rigorously justified 
it does represent the "state of the art" . 

Having developed the general framework, it is now specialized to the case of strong 
dissipation due to the bath. This assumption greatly simplifies the analysis as it means that 
(a) the dissipative stresses can be neglected since they will act more slowly than the friction 
due to the bath and (b) the momentum current will always be small so that the quadratic 
convective term can be neglected. The result is 

^ + iv.j(r)=C (11) 
at m 

dt + 9 (r) V Ji$) = " 7j (r) + (r; *) 

The linearity allows us to eliminate the momentum giving 



i + aw = Mv . u )y %*m _ v , (r;t) (12) 



7 dt 2 dt 7m \ 5p (r) / y 7m 



Scaling time by I/7 shows that the second-order time derivative is of higher order in I/7 
than the other terms so that it can be neglected in the strong damping approximation leaving 



dp(r) k B T„ ( , ,„5/3F[p)\ „ 2k B T .... . 



dt 7m \ 5p (r) J y 7m 

(Note that this result can be obtained directly from the full hydrodynamic equations by 
introducing a scaling of the time and momentum current of t = ^t* and j = ^ppj*, with p a 
typical density and I a molecular length scale, and noting that the dissipative stress tensor is 
of at least first order in the momentum current.) Although the coefficient of the noise term 
is state-dependent, it is nevertheless the case that the Ito and Stratonovich interpretations 
of this stochastic differential equation are the same as shown in Appendix |A] This model is 
well-known in the literature as the Dean-Kawasaki model^i^ and it has been widely used, 
e.g., in the study of the glass transition. Chavanis^ discusses various approaches to its 
interpretation. The derivation sketched here emphasizes the relation between this model 
and the more general theory of fluctuating hydrodynamics because the latter can be used, 
in conjunction with the approach described below, to address related problems to which 
Brownian dynamics is not applicable such as nucleation in simple fluids, multicomponent 
reacting fluids, etc.. However, for the specific problem of nucleation of over-damped colloids 
it is worth noting that the Dean-Kawasaki model can be derived by more rigorous methods 
as discussed e.g. by Kawasaki anc l Chavanis^. 

Without the noise term, this is often referred to as Dynamical Density Functional Theory. 
Another way of thinking about this is to say that the usual DDFT results from averaging 
this result over the noise. The question then arises as to whether the "free energy" entering 
here is the thermodynamic free energy or some other quantity (such as a coarse-grained free 
energy) that becomes the free energy after averaging over the noise. From the sketch given 
above, it would appear consistent to interpret the free energy functional in Eq.f JTB"]) as being 
the free energy functional for the one-component system of large molecules or colloids since 
it arises from a local-equilibrium approximation for the pressure which is in turn entirely 
due to interactions between the large-particles. 

Note that at low density one has^ 

0F[p) ~ y"(p(r)lnp(r) -p{r))dr (14) 



so that the dynamical equation becomes 



dp(r) k B T 2 2£; B T 

777 Vp r -V-eW p(r)€(r;t) 15 

at 7m y 7m 

which is simple diffusion with noise. One can therefore identify D = as the diffusion 
constant for the large particles in the bath so that it, and therefore the coefficient 7, can be 
determined from the Einstein- Stokes equation based solely on the size of the particles and 
the properties of the bath. 



B. Spherical symmetry 

In principle, this model could be studied numerically. Here, however, the goal is to make 
contact with the standard approach to nucleation and one of the key assumptions usually 
made is that of spherical symmetry of the clusters. One cannot, however, simply assume a 
spherically symmetric density profile in Eq. (j!3p as the noise term breaks spherical symmetry. 
One solution would be to average Eq.( fl"3]) over a spherical shell but it proves more convenient 
and natural to integrate over a spherical volume. 

The total (or cumulative) mass inside a spherical shell of radius r, m(r), is simply the 
integral of the density 



m (r) = / p (r') dr'. (16) 

J r'<r 

Integrating the dynamical equation over a spherical shell gives an evolution equation for the 
cumulative mass 

^ = D l V. [p (,0 Vffl) dr' - e [ V- ^Dp~V)i (r'; t) dr' (17) 

= D [ (p(r) V S/3 /}d ) rdr-e [ V-^2Dp (r% (r' ; t) dr' 
J V op{v) J J r , <r 

where Gauss' theorem has been used to convert the first integral on the right to a surface 

integral. It is straightforward to show that the autocorrelation of the noise term is 



6 (r - r') S(t- t') 2DeHnr 2 (^-^ J 5 (r" - r) p (r") dr"^j 
where the last term is the density averaged over a spherical shell of radius r, 



(18) 



<P]R>S ^4^ I S (r-R)p(*)dr (19) 
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Hence, the noise term can be replaced by e^8nDr 2 < p;r >s£(r; t) where < £(r; t)£(r'; t') >= 
5 (r — r') 8 (t — t'). At this point, since there is no more coupling to the non-radial noise 
terms, one can consistently seek a solution in terms of spherically symmetric density distri- 
butions, p (r), thus giving 

dm (r) 2 d 5/3F [p] 

-DAnr p{r) — 



-ey/8Trr 2 Dp(r)£(r;t). (20) 

p(r) 



dt dr 5p (r 

It is again easy to verify the equivalence of the Ito and Stratonovich interpretations of this 
equation (see Appendix |A|) . 

This result can be made somewhat more transparent by noting that for a spherically 
symmetric function p (r) = p (r) one has 

P (0 = J 5(r -r')p (r') dr' (21) 

so that for any functional F of p (r) 

8F [p] _ r SF [p] 5p (rO ^ = SF\p) 1 



<5p (r) y <5p (r') <5p (r) <5p (r) 4-7rr 2 

giving 



9m fr) 



Next, using the functional chain rule gives 

d ( 1 5pF[p}\ d 1 Z" 00 tfm(r') 5(5F [p] 



dr \47rr 2 Sp(r) J dr 4irr 2 J 5p(r) 5m (r') 

(r - r) Airr ; , 1 J dr 



dr' (24) 



dr Airr 2 J Q 6m (r') 

5m (r) ' 

where Q(x) is the Heaviside step function, so that the SDE can be written as 



«9m(r) = _ D d_np_ _ / 2£) 9m(r) 
ot ar om (r) V ar 

This shows very clearly that the spherically symmetric dynamics is gradient-driven with 
kinetic coefficients D drr !^^ and obeys a fluctuation-dissipation relation. It also emphasizes 
the central role played by the mass in this approach to nucleation - as opposed to density 
or any order parameter - which is a consequence of the fact that the theory is founded in 
hydrodynamics where conservation of mass is always true. 
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C. Nucleation pathway 



The stochastic evolution equation, Eq. (125j) . could be used as the basis of numerical sim- 
ulations. However, the goal here is to provide a description of the nucleation pathway and 
to make connection with more phenomenological descriptions of the nucleation processes. 
It is therefore necessary to characterize the typical or expected pathway. One natural way 
to do so is to determine the most likely path (MLP). When the noise is weak, it can be 
expected that typical paths will be close to the MLP so that the latter can be viewed as 
characterizing the process. If the noise is strong, the status of the MLP is less clear and it 
may only be one of many alternative paths that could be explored by the system. 

1. The path probability 

The question of the determination of the most likely path for a stochastic process dates 
back to Onsager and Machlup^i. For Gaussian processes with constant diffusion matrix the 
result is relatively easy to derive whereas for the most general cases there are several sub- 
tleties and the most general result seems to have first been derived by Graham^. Here, and 
in the following sub-subsection the theory of the MLP is reviewed and then its application 
to the problem of nucleation is discussed. 

Consider the general stochastic dynamics given by 

(IT • 

^ = b i (x(t)) + eQ ij (x(t))Z j {t) (26) 

where the noise terms, £ • (t), are Gaussian, white and uncorrelated. Note that here, and 
in the following, the Einstein summation convention is used. Interpreted within the Ito 
calculus, Graham^ shows that the probability to observe a given continuous path is given 
by 

P (x (t)) ~ exp ^- jT N C (x (t) , x (t)) dtj (27) 

where the Lagrangian is 

C (x, x) = \e~* ( Xi - C{ (x)) Df (x) ( Xl - Q (x)) + * det (Q (x)) + ^ 

(28) 
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with 



Da = QuQji (29) 
c ^_I e 2 det(Q) 9 



2 dxj det(Q)' 

the Ricci curvature scalar formed using D^ 1 (x (£)) as a metric in x— space is 

R = D ab D cd l f + d ^ d ^ d ^ 



'2 \dx c dxd dx a dxt, dx a dxd dx b dx c 

+D ab D cd D; f l (r: b r f cd - q 6 rt) 

and the Christoffel symbols of the second kind are 

The probability of observing a path is simply the probability of observing the required noise 
needed to generate the path plus additional terms related to the change in variable from the 
noise to the path. The first term on the right in the Lagrangian is due to the fact that the 
probability to observe a sequence of noises is by definition P (£ (t)) = exp f-± j t [ N f (t) dt) 
together with a replacement of (t) based on rearranging the stochastic differential equation. 
The second and third terms are due to the Jacobean in transforming variables from £ (t) 
to x(t). Note that the first term is of order e -2 , the second of order e° and the third of 
order e 2 so that in the weak noise limit, only the first term contributes. In this case, another 
interpretation of the expression for the path probability is that it is the Freidlin- Wetzel 
action functional from the theory of large deviations' 1 ^. The fact that these results can be 
written in the fully covariant form given above is, as noted by Graham', necessary since the 
path probability must be independent of the choice of variables used to express it. Thus, 
while the original stochastic model has no explicit geometrical interpretation, the focus of 
paths, which are geometrical objects, induces a natural geometrical perspective in which the 
(inverse of the) covariance of the noise plays the role of a metric in the space of stochastic 
variables. 

Given fixed endpoints, the MLP is determined by maximizing the probability with respect 
to variations in the path resulting in the Euler-Lagrange equations 

d d „ / , , ■ , a d 



-C (x (t) , x (i)) - — C (x (t) , x (*)) = (31) 



dtQ x . \ J dx 
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or, after some manipulations, 



d 2 1 dx c dxi 
dV Xr+ 2 lc dt dt 



D r 



dxi 



D r 



d_ 

dxi 



dx c 

Cm 



dt 



(32) 



r^ +e 1 det(g) ^^(Q) 



R 
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This can be interpreted as describing a particle with curvilinear coordinates x moving in 
response to a force given by the gradient of a potential under the metric D _1 . (Note that 
here and below, it is assumed that one selects a solution that actually corresponds to a 
maximization of the probability, as opposed to a minimum or saddle point extremization.) 

Equation (l32|) was derived under the condition that the end points of the path are fixed, 
so it defines a two-point boundary value problem with fixed initial point, x (0) = x , and 
end point, x(tf) = x/ and the time, tf, is determined by minimizing the action . In some 
cases, the endpoints are not both known a priori and one wishes instead to determine the 
MLP from any state satisfying some set of constraints to any other satisfying a different 
constraints. (For example, in the nucleation problem, this might be the path connecting a 
cluster of total excess mass Mi to one with total excess mass M 2 .) If these constraints involve 
only the coordinates and are represented in general as K a (x(£ )) = and K a (x(t/)) = 
then the boundary conditions become 



i, (0) = Ci (x (0)) + £ X^DJ (x (0)) (^M.) (33) 

a V ox i / x (0) 



x(0) 

dLg (x) 
dxi 



where the A's are Lagrange multipliers which must be determined. To illustrate, suppose the 
desired path started at point x but that the final point was determined as having a fixed 
mass, M(x.(tf)) = Mf. If there were n parameters, there would be n constraint equations, 
Eq.f l33p at the end point Xf together with the condition on the final mass. These n + 1 
constraints would take the place of the n conditions otherwise determined by specifying the 
final coordinate as well as providing one additional constraint for the Lagrange multiplier. 
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2. The most likely path when there is a fluctuation- dissipation relation 



In the limit of weak noise, i.e. for e — ?• 0, it can be shown (see Appendix |B]) that for for 
the case of a generalized potential driven dynamics of the form 

fc(x) = -L^)%^ (34) 

that obeys a fluctuation-dissipation relation, 

Lij(x) oc Aj(x) (35) 

the MLP between two metastable states, i.e. points satisfying foj(x) = 0, must pass through a 
saddle point and follow the deterministic path connecting the saddle point to the metastable 
states, 

-dt = ±L - (X) ^- (36) 
It is easily verified by substitution that this is an exact solution to the general expression, 
Eq.( )32|) . when Eqs. (I34H35P hold and the initial and final states are metastable points. It can 
be interpreted as saying that, starting in the metastable state, one follows the time-reversed 
dynamics against the free-energy gradient up to the critical point; beyond the critical point, 
one follows the forward-time dynamics down to the stable state. Note that this is only a 
mathematical prescription that serves to determine the MLP: there is no suggestion that 
the highly dissipative physical system is in any sense time-reversal invariant. Similarly, the 
MLP is not the same as the deterministic (i.e. e — > 0) limit of the stochastic process: it is a 
separate concept that concerns the most likely path to move both up and down free-energy 
barriers whereas the deterministic path can only describe movement down a free-energy 
gradient and is a priori incapable of describing movement up an energy barrier. While it 
so happens that for a part of the path between two points that involves moving down a 
free-energy gradient, the MLP and the deterministic limit coincide, this is not a general 
identification of the two: thus, while the deterministic path is indeed given by Eq.( l36p with 
the minus sign, the more general - and most certainly fluctuation-driven - MLP is given by 
Eq.( )36|) with the appropriate sign for the part of the path being determined. 

Recognizing that Eq.( l36|) is equivalent to steepest descent in curvilinear coordinates with 
metric = L^ 1 means that the MLP can be determined using the standard gradient- descent 
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algorithms^. This involves first locating the saddle point, x s , using, e.g. eigenvector- 
following techniques^, and then perturbing slightly in the direction of the of the unstable 
eigenvector as determined from the generalized eigenvalue problem 

HijVj = XL^Vj (37) 

where the Hessian at the saddle point is 

d 2 v(x) 



(38) 

In principle, the perturbation should be infinitesimal but in practice, some small, finite 
perturbation must be used with the approximation becoming exact as the size of the pertur- 
bation goes to zero. From this initial point, the steepest-descent path is determined using 
Eq.( l3T)|) with the positive sign (i.e. the forward-time, deterministic dynamics). 

3. Application to nucleation 

For the nucleation problem, instead of a collection of variables, Xi, as considered above, 
one has the density p (r) or the cumulative mass function, m(r), which are both fields. 
Generally, stochastic differential equations involving fields are understood to be defined by 
some sort of spatial discretization scheme so that the field, say m (r), is actually a collection 
of values, m (r») = rrtj, defined on a lattice thus allowing the use of the results developed for 
discrete variables. 

With this understanding, it is easy to verify that the general model, Eq. flIB"]) . does indeed 
posses a fluctuation-dissipation relation between the deterministic and fluctuating parts of 
the dynamics so that one can immediately infer that the MLP between the initial and final 
phases will pass through a saddle point determined from 

5p(r) 

where \i is the chemical potential of the mother phase and will be determined by moving 
away from this state in the forward and backward directions either by following the time- 
reversed dynamics from the metastable state to the critical state and then the forward-time 
dynamics thereafter, 

1(3 



or, equivalently, by perturbing slightly away from the critical cluster in the direction (or 
anti-direction) of the most unstable eigenvector and following the forward-time dynamics 
to the respective minima. This is the first significant result of this analysis and represents 
a generalization of the usual Dynamic Density Functional Theory to barrier crossing prob- 
lems. Lacking this result, applications of DDFT have until now had to rely on the artificial 
introduction of noise via the initial conditions of the density in order to induce nucleation 



(for a recent example, see Ref. l36l ) . Such a procedure is in fact inconsistent in so far as the 
density appearing in DDFT is already a noise- and initial-condition-averaged quantity and 
so, for example, for mother phases which are fluids it should simply be a constant. 

In the case of spherical symmetry, the dynamical model given in Eq. (l25j) is obviously of 
the same form as Eq. fl26|) and obeys the fluctuation-dissipation constraints of Eqs. fl34ti35|) . 
The MLP in the limit of weak noise therefore passes through the saddle point determined 
from 



Q _ 50F [p] _ d 5f3F [p] 



5m (r) dr 5p (r) 



(41) 



p(r) 



and is given by 



dm(r) _^^dm(r) 5/3F [p] _^<9m(r) d 5/3 F[p] 



dt dr 5m (r) dr dr 5p (r) 



(42) 



p(r) 



which can be given a geometric interpretation as steepest descent on the free energy surface 
in mass-space with a metric of ( ~\ . Note that this is exactly equivalent to Eqs. (l39H40]) 
with the assumption of spherical symmetry, thus showing the equivalence of the approaches. 
It also allows one to define the distance between two mass distributions, m (r) and m\ (r), 



as 




dm(r,t) ( dm(r,t)\ 1 dm(r,t) , , / . 

s mo,?™-! =mmmm / \ drat (43) 

t paths J y Jo dt \ dr J dt v ; 

where the right hand side is minimized over all " acceptable" paths connecting the two mass 
distributions with m(r,0) = m (r) and m(r,T) = mi(r). An "acceptable" path is one 
obeying the basic constraints on the cumulative mass distribution that it be monotonically 
increasing as a function of r and that m(0) = 0. These properties mean that the position r 
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may be replaced by m(r) as the variable of integration yielding 




dr (m t*) \ 2 
min min / \l I I — - — ) dmdt 



T paths Jo V Jo \ Ot 



'0 



(m 1 (x) — m 1 1 (x)) 2 <ix 



where the last line follows from recognizing that the second line now defines a Euclidean 
metric. Since uiq 1 (x) is the location at which distribution mo (r) has value x and m^ 1 (x) 
is the analogous quantity for the second distribution, this expresses the distance between 
the distributions in terms of the sum over Euclidean distances between equal-mass points. 
This structure is entirely due to the presence of the factor drn Q^ ~ p (r) in the dynamical 
equations and, via the fluctuation-dissipation relation, in the noise. Without it, Eq.( T43]) 
would already be Euclidean and the distance function would simply depend on the integral 
over (m (r) — mi (V)) 2 . 

The fact that the exact distance function can be easily calculated means that powerful 
approximate methods for determining the steepest descent path, such as the String method^ 
and the Nudged Elastic BancP^, can be used to determine the steepest descent paths 
defined by Eq. (142 j) which is to say, the MLP. Another implication concerns the stochastic 
evolution equations themselves. The fact that the induced distance function is Euclidean in 
the variable r(m), immediately implies that if the SDE is written in terms of this variable, 
the noise will be white (i.e. the coefficient of the noise variable £(r; t) will be constant rather 
than state-dependent). That this is the case is easily verified by first observing that 

dr (m) dr dm f dm\~ l dm 
dt dm dt \ dr J dt 

Using the evolution equation, Eq. (l2"5]) . gives 



dr (m) f dm\ 1 dm{r) 5(3F[p] / / dm(r 



-i 



st I,* J -^r = ~ D 1^-i 2D {r^r ] f(r;i) (46) 

where the term must be expressed in terms of r(m). At the moment, it does not seem 
as if the goal of getting white noise is going to be achieved. However, note that one must 
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express the noise in terms of m as the dependent variable rather than r. Since 



(47) 



5(r (m) - r (m')) 5 (t - t') 



dr (m) 
dm 



8{m- rri) 5(t- t') 



it is clear that one can make the substitution £ (r; t) 



dr(m) 
dm 



£ (m; t) giving 



dr (m) 



S(3F[p] 



5m (r) 



e\f2D£ (m; t) 



(48) 



r(m) 



as expected. It is not clear whether this formulation is of any practical significance. 



III. ORDER PARAMETERS 



So far, the central quantity governing the description of nucleation has been the cumula- 
tive mass, m(r; t). However, the classical description of nucleation is typically formulated in 
terms of a small number, often just one, order parameter. For example, in Classical Nucle- 
ation Theory clusters of the new phase are assumed to be spherical, with a radius R, to have 
the bulk new phase in the interior of the cluster and to have an interface of vanishing width. 
The only quantity that can vary is therefore the radius and this is the order parameter in 
terms of which the dynamics are formulated. 

The point of view adopted here in order to make a connection with the concept of order 
parameters is that the parameters must determine an approximation to the cumulative mass 
distribution or, equivalently, to the density. For example, in the CNT picture described 
above, a cluster has interior density p which is just the bulk density of the new phase; 
exterior density which is the density of the metastable phase, a radius R and an interface 
of negligible width. Hence, it implicitly specifies a density model 

p ( r; t ) ■— ► p (r; R(t)) = p Q (R(t) - r) + Poo 6 (r - R(t)) (49) 

and a cumulative mass of 

4:71 47T 

m ( r; t )^m (r; R(t)) = —r 3 p o e (R(t) - r) + — (R 3 (t)p Q + (r 3 - R 3 (t)) Poo ) 9 (r - R(t)) 

(50) 
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In general, the approximating function can depend on many parameters so that the general 
case is p(r;t) i — > p(r;x(t)) for a set of time-dependent parameters x(i). The goal in this 
Section is to determine how the time-dependence of the restricted description in terms of the 
order parameters can be developed from the general theory. The reason for investigating this 
question is that it has been shown that one can derive the CNT description of the critical 
cluster as the lowest-order result in a systematic expansion of the DFT free energy based 
on these types of parameterizations of the densityi^^ so that this gives a path for making 
contact between the general theory and CNT. Finally, it is also relevant to note that any 
practical scheme for the integration of Eq. (|2"5|) will involve some sort of representation of the 
field in terms of a finite number of parameters, so the results obtained here will be of equal 
use in numerical calculations as discussed below. 

In the following, two different approaches to the use of order parameters will be described. 
The first is solely concerned with the determination of the most likely path for nucleation 
and involves relatively few assumptions. The second is aimed at determining approximate 
stochastic equations for the order parameters which, in turn, can be used to determine 
the MLP. The approaches have different strengths and weaknesses: the former involves 
fewer assumptions but is more restricted in its use while the latter gives broader contact 
with earlier work but is somewhat heuristic. Both methods become exact in the limit of 
a complete set of order parameters and of course agree with one another in that limit, as 
demonstrated below. 



A. Direct determination of pathway using order parameters 

In the weak noise limit, the probability for a path based on the spherically-symmetric 
dynamics, Eq.( 125|) . is 

p ~ exp (-skf M ) (51) 

with Lagrangian 

„ 1 f°° 1 (dm(r) „ a , , d 1 5f3F[p}\ 2 , 

It is shown in Appendix Q that if the density is parametrized as p(r;t) = p(r;x(i)), then 
the weak-noise Lagrangian becomes 

1 dx a dx b dx a d(3tt[p} 

c = 2 9 « b (x) -ww + D ir^r + v (x) (53) 
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where Q[p] = F[p] — pp and 

V(x) = -D 2 
v ; 2 

and where the metric g a b (x) is 

9ab 



r 2 p (r) 



d 1 SfiF [p] 
dr r 2 5p (r) 

dm(r) dm(r 



dr 



dr 



o Airr 2 p(r) dx a dxb 
Maximizing the path probability gives the Euler-Lagrange equations 

dgu dg i:j dg t / 



9n 



d xi , 1 dx\ dxj 



dt 2 



+ - 



or 



2 dt dt 



dxj dxi dxi 



dxi V X 



da;; fix. 



dt 2 1 ' Zj (it dt 

where the Christoffel symbol of the second kind is 



1 



kl 



dg m k dg m i 



dg 



ki 



(54) 



(55) 



(56) 



(57) 



(58) 



dx t dx k dx m / 

Equation f )57|) is recognized as the equation of motion of a particle expressed in curvilinear 
coordinates with a force derived from the potential V(x). If the end points are determined 
by the constraints, 



J,(x(0))=0, Z = l, 
^(x(T))=0, / = 1, 

then the boundary conditions are 



n 



n 



(59) 



dxi 


+ 

t=o 




<9/3F [p] 




~dt 


dxi 


t=o 


dxi 


+ 

t=T 




d(3F [p] 




~dt 


dxi 


t=T 



At/ (x; 



9J r 



(60) 



t=o 



where /x r and A r are Lagrange multipliers introduced in the course of minimizing the action 
under the constraints. Note that if, e.g., the initial point is fixed, then J r = x r (T) — xo r 
and the constraint just serves to evaluate the (uninteresting) Lagrange multiplier: in this 
case the constraint equations play no role. On the other hand, if a nontrivial constraint is 
applied then instead of fixing the value of x r (T) a condition is imposed on the velocities at 
time T. For example, for nucleation one might know the initial state, x (0), and might want 
to know the most likely final state with a given excess mass, 



M = 4?r 



(p (r;x(T)) - Poo )r 2 dr. 



(61) 
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Then, the boundary conditions at the end point would be 



dxi 
~~dt 



D _ u ,df3F[p] 



t _rp <±7I UX\ 

B. An order-parameter dynamics 



t=T 



The spatial and temporal components of the density can be separated by writing it in some 
parametrized form as p(r;t) = p(r;x(i)) where x(t) is a collection of scalar parameters. 
One common method of doing so would be to expand the density in a complete set of basis 
functions. Alternatively, the parameters could represent a discretization of the density, 
x n = p (nA) for some small distance A. As these examples indicate, there will be in general 
an arbitrary number of parameters. Their time evolution can be developed by noting that 
the m(r) will also be a function of the parameters so that from Eq. (1201) one can write 

dm(r;x(t))dxi 2 d 5f3F[p] 

- UAnr p (r; x {tjj 



-e^8nr 2 Dp (r;x(t))f (r;t). (63) 

p(r;x(t)) 



dxi dt ' dr 5p (r 

Suppose that this is multiplied by some function Wj (r; x (t)) and integrated giving 

80 (x) t = D f w ' (r; x (t)) w " (r; x (t,) I ( w) ^ dr (64) 



Wj (r; x (i)) e^8nr 2 Dp (r;x(t))f (r; t) dr. 



where 



0y (x) = jH Wi (r; x (*)) ^(^(*)) dr (65) 
The autocorrelation of the noise term is easily calculated giving 

/■OO 

e 2 2D / 47rr 2 W i (r;x(t))W 3 - (r; x(t)) p (r; x (t)) dr (66) 
Jo 

While these equations are exact regardless of the choice of Wi(r;x(i)), it will soon be 
apparent the the result is particularly simple if the diffusion matrix and the matrix (x) 
are required to be proportional which is only possible if, modulo a multiplicative constant, 

W - (v x (t)) = - dm(r;x(t)) 

W>V>XW) 47rr 2 p(r;x(t)) dx . ^0 

so that T>ij (x) = 2Dgij (x) with 

r 1 0m(r;x(t))0m(r;x(t)) 

9vW ~J 4vrr2p(r;x(t)) dx, dx, ^ { ™ } 
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i.e., the same metric as found above (see Eg. (1551) ). Another result of this choice is that the 
thermodynamic driving force becomes 

d 8f3F[p] 



W,(r;x(t))47rr^(r;x(t)) 6 
<9m(r;x(t)) d 5(3F[p] 



dr 



(69) 



p(r;x(t)) 



dxi dr Sp(r 
dm(r;x(t)) 5^F[p] 



dr 



p(r;x(t)) 






oo 


p(r;x(i))_ 


.-I. 




oo 


p(r;x(t))_ 


.-/ 



«9 2 m(r;x(i)) 5/3F [p] 



dxidr 8p(r 
dp(r;x(t)) 6f3F[p] 



dr 



dxi 



dp (r) 



p(r;x(t)) 

dr 



p(r;x(t)) 



For the problem of nucleation, the first term gives no contribution at r = while for very 
large r, one expects the system to have the bulk properties giving 



lim 

r->oo 



dm(r;x(t)) 6pF\p] 



p(r;x(t)) 



dN (x (t)) 
dxi 



-fi 



(70) 



dxi Sp (r) 

where is the total number of particles and [i is the chemical potential in the bulk. Using 
the functional chain rule, the second term is recognized as 



0p(r;x(t)) 5(3F[p} 



dr 



p(r;x(i)) 



d/3F (x) 
dxi 



(71) 



dxi 5p (r) 

where F (x) = F [f\. (Note that this is only exact if the parametrization is complete.) 
Combining these results gives the stochastic model 

dm (r; x (t)) 



9u 



-D 




2D 



-£ (r; t) dr. 



(72) 



dt dxi J \l Anr 2 p (r; x (£)) 

F — //AT. Ideally, one would like to replace the noise by a simpler form 



where again Vt 

giving the same autocorrelation function however there is a complication. The use of the 
usual chain rule for derivatives means that the resulting stochastic differential equation must 
be understood in the Stratonovich interpretation^ and, unfortunately, in the present case 
the resulting equation is not Ito-Stratonovich equivalent. It is therefore not the case that 
one can substitute one noise term for another with the same autocorrelation matrix as this 
only holds in the Ito form. As described in Appendix [Fj the spurious drift arising from 
Ito-Stratonovich inequivalence gives rise to an additional contribution to the deterministic 
part of the equation with the final result 



% = (x) + 2De>A, t (x) - eV2D q ^ (x) ^ (t) 



dt 



dxi 



(73) 
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with gij = qaqji (note that g is positive semi-definite so this decomposition is always pos- 
sible) and what will be termed the "anomalous force" (since it does not arise from the 
thermodynamic driving force) is 

, dq~} (x) 1 . dq^t (x) 
A i (x) = fc 1 (x) - - 2 gu l (x) -^^9^ (x) (74) 

+- (a" 1 U) a' 1 U) - a" 1 fx) a" 1 U)) f° - ^x) dm (r; x) dm (r; x) 

The presence of the anomalous force means that there is no exact fluctuation-dissipation 
relation. However, since it has its origin in the noise, and is therefore of order e 2 , it will not 
affect the weak-noise limit and the MLP will be determined as usual from the forward- and 
backward-time dynamics starting at the critical point. Again, it is important to remember 
that in the strong-noise regime Eq. ( I73p must be understood in the Stratonovich interpre- 
tation : the equivalent Ito equation will have a modified form of the anomalous force (see 
Appendix [F] for details). 

The advantages of the choice for Fj (r; x (t)) are now evident. First, the particular choice 
used here results in the free-energy term having the form of a simple gradient of the free en- 
ergy due to the ability to integrate by parts in Eq.()69l). Second, the fact that the fluctuation- 
dissipation relation is preserved means that in the weak-noise limit one recovers the classical 
behavior with the MLP passing through the critical cluster. 

Equation f!73|) is exact when the parameters are determined from a complete represen- 
tation of the density. For a finite collection of parameters, it seems reasonable to continue 
to use Eq.( 173i) as an approximation and numerical evidence in support of this will be given 
below. Appendix [E] gives an alternative derivation of this model showing that it is reason- 
able even for an incomplete (i.e. approximate) representation of the density provided that 
the difference between p (r; x (t)) and the actual density p (r; t) is in some sense small. 



C. Most likely path 

A heuristic stochastic dynamics for a set of order parameters is given in Eq.( |73l) above. 
In the weak noise limit, in which the anomalous force can be neglected, this dynamics is 
gradient-driven and satisfies a fluctuation-dissipation relation (FDR) so that it immediately 
follows that the MLP connecting the two metastable phases passes through the saddle point 
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defined by 

d/3Q , . 

0= ik < 75) 

and is determined by steepest-descent according to 

dxi _! <9/3fi 

dT = ^ (x) ^7 (76) 

Following the same arguments given above, the general MLP connecting any two states is 

determined from the Euler-Lagrange equations 

d 2 xj 1 fdg is dg ir dg rs \ dx r dx s d ( D 2 dj3Vt _ l df3Vl \ 
9ij dt 2 2 \dx r dx s dx t J dt dt ~ dxi V 2 d Xj 9jl dx { J 1 j 

It is easy to show by direct substitution that Eq. (17(51) is a solution to Eq. ( 1771) . 

Two expressions for the evolution of order parameters have been derived: Eq. fl76|) . which 
follows from the general equations for the MLP, Eq.( l77[) . for the special case that the end- 
points are metastable states, and Eq.f l56p . which follows from the minimization of the path 
probability after assuming a particular test function. In essence, the former results from 
introducing the approximate density after determining the equations for the MLP while the 
latter results from introducing the approximate density before determining the MLP. Each 
formulation has its advantages: the gradient formulation shares certain properties with the 
exact MLP such as the role of the critical cluster whereas it is not even clear a priori that 
Eq.( )56|) will give a path passing through the critical point. On the other hand, Eq. ([56|) is 
exact, given the assumed form for the density, whereas the gradient equations are derived via 
a series of manipulations which are not unique: one could derive any number of equations 
of motion for the parameters from Eq.( l63l) by similar manipulations to those used above. It 
is therefore of some interest to examine the connection between these two formulations. 

Comparing Eq. (1571) and Eq. (177j) it is clear that the only difference is in the source terms 
and that equivalence would demand that 

°° r 2 p (r) (±L 6 1EM) 2 dr - -i^ (78) 

rp[r) \drr 2 5p(r) J dr ~ d Xj ^ 1 dx L [ ™ } 

Now, from the derivations above, it is already known that 

mi = _r™ dm(r) d 1 6PF[p] dr 
dxj J dxj dr r 2 5p (r) 

so 



P t P = r dr , r drr i / r) (9 1 S(3F[p] \ / d 1 6PF[ P ] \ f_ 21 dm(r) ^ dmir') 
dxj jl dxi J Jq \drr 2 5p{r) J \dr' r' 2 5p{r') J \ dxj jl dxi 
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which is the same as Eg. (ITS]) provided that 

^'^=^W«M. (80) 

As demonstrated in the Appendix O, this is a completeness relation and is sufficient to 
give equality between the two descriptions of the MLP. It is satisfied when the parameters 
result from an expansion of the density in a complete set of basis functions. The use of the 
order-parameter dynamics for a finite collection of order parameters is therefore heuristic. 
However, as shown below, only a small number of parameters is necessary for the heuristic 
dynamics to show convergence to the exact result. 



IV. APPLICATIONS 

The classic problem of the nucleation of a liquid droplet in a supersaturated vapor will 
be used to illustrate the theory described above. Calculations were performed using the 
squared-gradient free energy model, 

F[ P ] = f (f(p(r)) + ^(Vp(r)) 2 ) dv (81) 

where f(p) is the bulk free energy per unit volume. The calculations reported here were 
performed for the ten Wolde-Frenkel model potential for globular proteins^ which consists 
of a hard core of radius a and an effective pair potential outside the core, 



<r>°) = -A\ " « 77T^ (82) 




where e is the energy scale of the potential and the parameter a controls its range. For large 
values of a, such as a = 50 used here, this gives a phase diagram typical of proteins including 
a metastable high-concentration ("liquid") phase^ 1 ^. The results presented here pertain to 
the transition from the low-concentration ("gas") phase to a high-concentration ("liquid") 
phase which is known to sometimes play a role as the first step in the processes of protein 



crystallization^ 1 ^. The bulk Helmholtz free-energy function, / (p 
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was approximated using 



and 
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The coefficient of 



first-order perturbation theory*^ as described in detail in Refs. 
the squared-gradient term, K, was evaluated using a recently-derived approximation based 
on the potentiate. 
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The most direct approach towards the determination of the nucleation pathway based on 
the field equations, Eq. fj4Tj) and Eq.( l4"2]) . would be to discretize them by first discretizing the 
radial coordinate, r i — > = 5, for some small length 5, which then implies a discretization 
of the density, p (r) i — > Pi = p (i5). Since the density is only expected to vary quickly near 
the interface, the efficiency of this scheme can be significantly improved by allowing for a 
dynamic discretization whereby the density and position of the lattice points are allowed 
to vary. Here, this is done by introducing a particular model for the density based on a 
continuous, piecewise-linear approximation, 



p[r) 



Po + (Pi - Po) 



p , r < w 

™, w < r < w + w 1 



Pi + (P 2 - Pi) 



W[) — Wl 
W2 ' 



w + wi < r < wq + wi + w 2 



(83) 



Pc 



w + ... + Wtv-i < r 



The parameters that are allowed to vary freely and dynamically are the combinations p { , Wi 
(referred to as " links" ) while the density p^ is always fixed at that of the surrounding bulk 
(e.g. the vapor density for nucleation of liquid droplets). Note that the simplest model 
consisting of just two links, 



p{r) 



Po + (Pc 



p , r < wo 
- po) ^p, w <r<w + w 1 
Poo, wo + wx < r 



(84) 



gives the minimal description of a cluster in terms of a radius, say wq, an interfacial width, 
W\, and a central density, p . Similar parameterizations based on these three quantities 
are common in the literature and here I note in particular one introduced by da Gama 
and Evans^ and used recently by Ghosh and Ghosh^ to study both planar interfaces and 
spherical clusters, 



p{r) 



Po 



Po - Pc 



exp (a (r — R)) 



®{R-r) + 



Q(r-R) 
(85) 

where the radius is here denoted R. The parameters occurring in all of these approximations 
are "order parameters" in the sense discussed above and the most likely path of the density 
can be determined by means of the order-parameter equations, Eq. fl76|) . These are of course 
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only approximations but the piecewise-linear scheme becomes exact as the number of links 
is increased. 

As discussed above, once a parametrization is chosen, there are three steps involved in 
determining the MLP. The first is to determine the critical cluster, i.e. the saddle point 
of the free-energy surface in parameter space. This was done using standard eigenvector- 
following methods^ 1 ^ which result in a given set of parameters describing the approximation 
to the critical cluster, x c . At the saddle point, the Hessian matrix determined from the free 
energy has one negative eigenvalue and a corresponding eigenvector, v ncg , which specifies 
the unstable direction. The second and third parts of the calculation are the numerical 
integration of the order-parameter equations, Eq. (l76l) . with initial conditions consisting of 
a small perturbation away from the critical point in the unstable direction, x c ± ev neg where 
the parameter e is chosen so that the change in free energy is some specified amount. The 
result of one of the perturbations will be that the system falls back to the metastable state, 
i.e. the cluster will diminish, while the other will send the system towards the final, stable 
state which means the cluster will grow. 

A. Nucleation of droplets 

The critical cluster for temperature ksT/e = 0.375 and supersaturation S = p v /p vc = 
1.175, where p vc is the vapor density at coexistence, was determined for both the piecewise- 
linear profiles with varying numbers of links and the exponential profile. The "exact" critical 
cluster was determined by solving 



numerically using a relaxation technique^ where the chemical potential is simply /i = /' (p v ). 
The excess free energy and excess mass of the critical cluster is given in Table 1 as calculated 
from several approximations. It demonstrates the gradual convergence of the properties of 
the critical cluster to the exact result as the number of links in the piecewise-linear profiles 
is increased. The exponential profile gives a relatively good description of the critical cluster 
based on only three parameters. Figure 1 shows the exact density profile of the critical cluster 
and the piecewise profile based on 2, 4 and 10 links. It is apparent that while the simple 
profile based on 2 links is quite crude, the 10-link profile is already a good approximation 
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TABLE I. Excess free energy, Af2, and mass (number of molecules, A A) for different density 
models. The label "PWM" refers to piecewise-linear profiles with M links, Eq. (|83p . 



Approximation /3Af2 AA" 

PW2 98.9 1471 

PW4 82.9 1250 

PW6 79.4 1211 

PW8 77.8 1186 

PW10 77.1 1175 

Exponential, Eq.jHSD 80.4 1246 

Exact, Eq.(j86D 75.8 1158 
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FIG. 1. The density profile of the critical cluster as determined by numerically finding a stationary 
solution of Eq. (|86p . the "exact" solution, and from the piecewise linear model with A links, Eq. (|83p . 
The inset shows the convergence of the central density towards the exact value as the number of 
links increases. For the model with 10 links the figure shows the boundaries of the links as circles. 
Notice how the density of the discretization is automatically adjusted to be highest in the region 
of most rapid variation of the profile. 

to the continuous profile. 

The evolution of the droplet was determined beginning with the profile for the critical 
nucleus. The Hessian for the free energy was evaluated and diagonalized giving in all cases 
a single negative eigenvalue. The profile was perturbed by adding a multiple of the cor- 

29 



< '■. ' 0.2 0.4 0.6 0.8 




2 



0, 





1 1 


1 


N/N 




p(0)/p c (0) 












AQ.IAQ.~~ 

C 






1 





250 500 750 



s 

FIG. 2. The variation with path distance, s as defined in Eq. (|87|) . of the excess number of molecules, 
N, total excess free energy, AO, and central density, p(0), relative to the values for the critical 
cluster for a nucleating droplet. The calculations were performed using the piecewise- linear model 
with 10 links (19 parameters). For reference, the values for the critical cluster are TV" = 1175, 
A/30 = 77.1 and p(0)<7 3 = 0.68. To aid in comparison, the vertical line picks out the same value 
of s as in Fig. ([3]). 

responding eigenvector with the coefficient chosen so as to give a change in free energy of 
lk-^T. The equations for the MLP, expressed in terms of order parameters, Eq. (!76|) . were 
then integrated using the Intel ODE solver library^. Depending on the sign of the coefficient 
used to perturb the state the system would either fall backwards towards the metastable 
state (i.e. decreasing cluster size) or forwards towards the stable state (increasing cluster 
size). The integration was stopped in both cases when the excess free energy reached lk-^T . 
For the forward (growing) direction, this was simply a convenient place to stop but in the 
backward direction, some such cutoff is necessary since the weak-noise assumption must 
break down. 

The independent variable in the integrations is time but care must be taken in interpreting 
the results. The "time" has the physical meaning of the time taken by the deterministic 
dynamics to drive the system in the desired direction. When moving forwards in time, this 
is therefore the physical time. However, when evaluating the MLP via the time-reversed 
dynamics, this "time" does not correspond to the physical time needed for fluctuations to 
drive the system up the free-energy gradient. Discussion as to what physical meaning can 
be attached to this time can be found in Bier et al^. In order to avoid any ambiguity, the 
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FIG. 3. The variation of the equimolar radius, R eq , and the total radius, Rtotai, as functions of 
distance along the path, s as defined in Eq. i|87|) . for a nucleating droplet. To aid in comparison, 
the vertical line picks out the same value of s as in Fig. ([2]). 

results here are displayed in terms of the dimensionless distance along the path denned as 



with s = being taken to correspond to the smallest cluster and the sign being chosen so 
that s increases monotonically moving from the smallest cluster to the critical cluster and 
on to post-critical clusters. 

Figure [2] shows the evolution of the central density, the excess mass and the excess free 
energy of the cluster as determined using the piecewise-linear parametrization with 10 links 
(19 parameters). At first, the central density remains very close to that of the vapor and 
the excess free energy is small - on the order of a few ksT. Surprisingly, even at the very 
beginning of the process the excess mass is finite and a significant fraction of the mass 
of the critical cluster. After some time during which the density and mass increase very 
slowly the growth enters a new regime in which the density increases rapidly to nearly 
the bulk liquid density. From that point onward, the evolution is unremarkable as the 
cluster grows to criticality and beyond. Corresponding behavior is seen in the size of the 
cluster which is characterized in Fig. |3]in two ways: the equimolar radius, R eq defined via 
y"-^eq(p(0) — Poo) = I (p( r ) ~~ Poo) dr, an d the total spatial extent, defined for the model given 
in Eq.( |84|) . as -Rtotai = X^o* w i- (For a general representation of the density, this might be 
characterized as the distance at which the density reaches some small threshold above the 




(87) 



31 



0.1 0.2 0.3 0.4 

^^links 



FIG. 4. The excess number of molecules in the cluster at the beginning of the nucleation process, 
N(0), as a function of the inverse of the number of links in the profile, iVi;^. The log-linear plot 
shows a nearly linear variation indicating that N(0) ~ Ae b / NliakB . Fitting all of the data to this 
functional form and extrapolating gives an asymptotic value of N(0) = 1331 and using only the 
three points with the largest value of A^n n k s gives N(0) = 1068. These values are consistent with 
the expectation that the limiting value should be the mass of the critical cluster, N c = 1175, shown 
as the square symbol. 

background.) Initially, the equimolar radius is nearly constant while the cluster has very 
large spatial extent and, counter-intuitively, the cluster shrinks . At the same point as the 
density begins its rapid increase, the equimolar radius also begins to increase even while the 
total spatial extent continues to diminish. Eventually, the latter reaches a minimum and 
the cluster grows according to both measures. 

These results are in stark contrast to the usual picture of droplet growth as assumed 
in the Becker-Doring picture (where droplets begin as dimers and grow by attachment and 
detachment of monomers) as well as the corresponding results from constrained DFT calcu- 
lations (where droplets begin as zero mass objects that grow monotonically) . In fact, how 
can a droplet "begin" with a finite mass? I shall attempt to answer this question from two 
different perspectives: first, that of the formalism and second from a physical point of view. 
In terms of the formalism, the MLP going from "nothing" to the critical droplet is the time- 
reversed evolution of the dynamics with no noise: in other words, it is the time-reversal of 
the deterministic evolution starting with a droplet slightly smaller than the critical droplet. 
What happens in such a case is obvious: the droplet evaporates by shedding mass to the 
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bulk. Since mass is conserved, the excess mass originally in the droplet cannot just disap- 
pear but must diffuse over ever larger volumes until it is lost in the (infinite) bulk. It is easy 
to understand, then, why the total mass is almost constant in the beginning of the process. 
Indeed, one would expect that the total mass should be the same as that of the critical 
cluster and should not change at all. That this is the actual result, if one could solve the 
model with the number of parameters tending to infinity, is indicated in Fig. [4] which shows 
the initial mass as a function of the number of links in the profile. It is clear that the mass 
increases with increasing refinement of the profile and that an extrapolation gives a mass 
close to that of the critical cluster. The fact that the mass is less than that of the critical 
cluster can be attributed to the difficulty of approximating the density distribution when 
the density becomes very dilute. Figure illustrates this as well, showing a large variation 
of the various parameterizations at short times with close convergence at long times. This 
includes the case of the exponential profile which is also not well suited to representing the 
diffusive distribution at short times. Note that according to this picture, one would expect 
the excess density at early times to behave purely diffusively, giving, e.g., 



(where "early times" means that t is large and negative). Then, the total mass would be 
conserved whereas the equimolar radius would be given by 



thus explaining the slow change in the equimolar radius. The fact that the result follows 
from such simple considerations suggests that, formally, it is very robust being independent 
of the assumed free energy model, potential or any other details. 

From a physical perspective, the explanation given so far is not very helpful. Instead, 
it is useful to divide the process into two parts as indicated in Figs. |2] and [3j During 
the first stage, matter comes together "from infinity" in the form of a diffuse structure. 
This can be understood as the formation of a long- wavelength density fluctuation. The 
presence of such fluctuations is expected and the fact that the structure formed only costs 
a few k B T in energy supports this expectation. What appears strange is the directed way 
that matter seems to come together but this is simply a result of two artifacts. The first 
is the assumption of spherical symmetry which means that the same thing must occur in 
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FIG. 5. The excess number of molecules, N, (left panel) and the total radius, R, (right panel) 
as functions of distance along the nucleation pathway using piecewise-linear models with N links, 
"PWN", and the exponential profile. The figures demonstrate that the short-time behavior is 
sensitive to how well the density profile is represented whereas all of the models converge at longer 
times. The exponential profile behaves very similarly to the PW2 profile, which has the same 
number of free parameters, but gives a better estimate of the properties of the critical cluster (see 



all directions at the same time; in reality, such a structure is expected to be messy and 
ramified. Second is the nature of the MLP. It is well known that the MLP generally follows 
the shortest possible path up the gradient, with no backtracking or deviation^. This gives 
it an appearance of being directed towards reaching a final state (i.e. the critical state) like 
a hiker deliberately climbing a hill. Any individual realization of the stochastic process will, 
in contrast, consist of moves in both directions - up and down the gradient - and will only, 
on average, eventually trace out the MLP (or something like it). A final artifact is that this 
is the MLP for a particular process: in other words, it represents a conditional probability 
in which the condition is that one does actually reach the critical state. A real system will 
make many abortive attempts to climb the potential gradient and only after many such 
attempts will it succeed. Here, in contrast, one studies the one successful attempt which 
again gives the appearance of being directed. 

So far, only the first part of the process - the formation of a density fluctuation - has 
been discussed. The second part of the process - marked by the rapid increase in density 
in the nascent cluster - can be interpreted as the occurrence of an actual nucleation event. 



Table U). 
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This part of the process essentially proceeds classically and can be visualized in terms of the 
Becker-Doring picture or of the typical DFT calculations. 

Putting these two pictures together gives a plausible and physically appealing interpre- 
tation. Density fluctuations of all types occur in the fluid and it makes intuitive sense that 
droplet nucleation should be more common in regions of higher density since the free energy 
barrier will be lower. Furthermore, the presence of a density fluctuation means that there 
is excess mass from which to build the cluster. Even in the classical picture (which does 
not assume any pre- nucleation enrichment) a cluster must somehow draw in material from 
the surrounding vapor but this part of the process is typically treated separately^, if at all. 
Considered in these terms, the classical process, with no enrichment, seems arguably more 
mysterious: a small, unstable cluster not only survives but manages to draw in material and, 
in the fluid, mass flows at such a rate that the cluster is always in contact with a vapor of 
exactly the bulk density. In the present picture, the needed mass is already locally present 
so that the formation of the cluster does not lead to the formation of a depletion zone or 
require the orchestrated flow of matter from far away. And finally, let us note that it is 
not the case that a fluctuation with exactly the right density appears in order to host the 
nucleation process: rather, fluctuations of all different sizes occur and if nucleation begins in 
one that is too small, then it may fail to reach the critical size. On the other hand, it is likely 
that the probability of a fluctuation decreases with its size so that the most likely scenario 
is of a successful nucleation event occurring in a fluctuation which is just large enough to 
supply the needed mass. This is the case that the MLP picks out. 

B. Nucleation of bubbles 

As a second application, the evolution of the reverse process - the formation of a low- 
concentration "bubble" within the high-concentration "liquid" - was determined. For these 
calculations, the temperature was the same as previously, ksT/e = 0.375 but the super- 
saturation was taken to be S = p v /p vc = 0.9875 so that the solution is under-saturated 
making the dense phase unstable towards nucleation of the less dense phase. The results 
are shown in Fig. ([6]). A bubble consists of a mass deficit relative to the bulk fluid and 
so, according to the deterministic dynamics, a sub-critical bubble is expected to vanish by 
gradually spreading the deficit over ever larger volumes. Therefore, based on the discussion 
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FIG. 6. The left panel shows the variation with path distance, s as defined in Eq. (|87p . of the excess 
number of molecules, N, the total excess free energy, Af2 and the central density, p(0). The right 
panel shows the evolution of the equimolar and total radii. The calculations were performed using 
the piecewise-linear model with 8 links (15 parameters). The excess number and free energies are 
scaled to the values for the critical cluster which are N = —546, Ail = 48.2&bT. 

of droplet nucleation, one would expect the process of bubble nucleation to begin with the 
formation of a long-wavelength, small amplitude deficit of density within which the bubble 
would nucleate. The results of the numerical calculation bear this out and all of the pre- 
ceding discussion of droplet nucleation can be carried over to the dual process of bubble 
nucleation. 




C. Order parameter dynamics 

As a final illustration of the formalism, the general expressions for the order-parameter 
dynamics, Eq.f J73|) . are specialized to give a CNT- level of description of the nucleation process 
for either droplet or bubble nucleation. By "CNT-level" it is meant that there will be only 
a single order parameter, the radius of the cluster R, and the calculations will be performed 
to lowest order in w/R where w is the width of the interface. The first step is the evaluation 
of the metric. This is most easily done for a concrete model and here the minimal piece-wise 
linear model with 2 links will be used. The central density is fixed at the bulk density of the 
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new phase, p , and the width w is held constant. Then, a straightforward calculation gives 



9rr 



dm (r) 



dr 



4tt(Po PcS R 3 {1 + 0{w/R)) 



(90) 



/o 47rr 2 p(r) \ dR J p ro 
where p^ is the density of the initial metastable phase. Using this, the SDE for the order 
parameter becomes 
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or, written in terms of the mass of a cluster, N = 4^i? 3 p , 



dN 
~dt 
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Recalling that this equation must be interpreted in the Stratonovich sense, it is equivalent 
to a Fokker-Planck equation for the probability of observing a cluster of size N at time t, 
p (N, t), of the formS 3 - 



dp (N, t) 
dt 



_d_ 

dN 



D 



4ttPoPc 



;R(N) 



d/3n 



p(N,t) 
(93) 



This can be compared to the classical result of Zeldovich (see Ref. Q, Eq. (9.27)) which is 
derived using the Becker-Doring picture, 



dp (N, t) d 
dt ~ dN 



p(N,t) 



(94) 



where the monomer attachment frequency, f(N), in the case of diffusion-limited kinetics is 
given by (see Reffl Eq. 10.18) 



f(N)=ri4nR(N) Dp c 



(95) 



and where all unknown details of the monomer attachment process are contained in the 
sticking coefficient, p. In the CNT limit of large N, the factor of 1/(6N) in Eq.( l93l) can be 
ignored and it then agrees with the Zeldovich result with the sticking coefficient identified 

as 



pI 



(PO - Poo)' 



(96) 



For the nucleation of droplets from vapor, one expects that p ^> p^ so that this gives rj ~ 1. 
This shows that by specializing the general theory to the CNT regime of large clusters and 
a single order parameter, one is able to recover all of the elements of the classical formalism. 



37 



Similarly, the free energy in the squared-gradient approximation is^2 

R 3 

n [p] - fioo = 4vr— (w - Woo) + 4rf 7 (97) 
where the excess surface free energy is 

and the density-averaged free energy per unit volume is 

1 f Po 

co = r / (u (x) - Woo) dx (99) 

[Po ~ Poo) J p x 

For a large interface, the width can be estimated by minimizing the free energy of the critical 
cluster giving 

» = J (100) 

where u cocx is the free energy per unit volume at coexistence^. All quantities here are 
therefore directly determined by the interaction potential. 
Substituting these elements into Eq. ( 1731) gives 



dR = _ D Poo (^0-^) ^-1 ( ! _ Re) _ D P°o(»°°-»°l R - H (t ) 

where 



(101) 



R c = — — — (102) 

is the usual CNT expression for the critical radius. Note that for a single order parameter, 
the anomalous force vanishes and so this model applies to the strong-noise regime. In 
principle, it can be used to determine a mean first-passage time and, hence, nucleation rate 
as will be discussed elsewhere. For now, it is interesting to make one further observation 
which is that for large, post-critical clusters, where the noise becomes unimportant, this 
gives the deterministic growth law 

d 4 « —D ^°° — — R" 1 (103) 

dt (P0 - Poo) 

implying that R ~ t 1 ^ 2 which is the well-known result for diffusion-limited cluster growth^ 1 ^. 
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V. CONCLUSIONS 



The goal of this work is a description of nucleation based on a dynamical, non-equilibrium 
description of fluctuations with the additional aim of making contact with the usual ideas 
of classical nucleation theory. The purpose is to develop a consistent and unique description 
of nucleation pathways, kinetics and post-critical cluster growth. Starting with fluctuating 
hydrodynamics this goal is achieved for the particular case of Brownian particles where the 
limit of strong dissipation provides several simplifications. The nucleation pathway is char- 
acterized as the the most likely path and a reduced description in terms of order parameters 
is developed. The order parameter description includes both a Langevin dynamics for the 
order parameters and a description of the most likely path. It is noted that in the weak noise 
limit, all levels of description lead to a unique description of the MLP that passes through 
the saddle point of the free energy functional and that is determined by steepest-descent on 
the free energy surface governed by a metric that is unambiguously specified. At this level, 
the goal of making contact with CNT is achieved with the role of the free energy emerging 
from the nonequilibrium description rather than being assumed a priori. Furthermore, as 
illustrated in the simple case of a single order parameter, one also has a Langevin descrip- 
tion of the process which can be used to develop a rate theory that recovers the classical 
results and is consistent with the expected post-critical growth law. Thus, the goal of giving 
a unified description of these three elements - rates, pathways and growth-laws - is also 
achieved. 

One of the main theoretical results concerns the determination of the nucleation pathway 
as characterized by the Most Likely Path. It is shown that, in the limit of weak noise, this 
is determined by the deterministic part of the dynamics moving either forward or backward 
in time away from the critical point. This result represents a generalization of Dynamic 
Density Functional Theory to problems involving barrier crossing. The same result holds 
when the dynamics are written in terms of order parameters. In all cases, the MLP can, 
alternatively, be viewed as being determined by steepest-descent on the free energy surface 
under a prescribed metric. It is noted that the strong-noise limit may also be handled but 
investigation of the effect of strong noise is deferred to a later date. In general, it is seen 
that one recovers the elements of CNT in the weak-noise limit which therefore plays a role 
analogous to the quasi-classical limit in quantum theory. 
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It is interesting to compare this result with previous approaches. As discussed in the 
Introduction, there have been two primary approaches to the description of nucleation path- 
ways using DFT. The first is the minimization of the free energy under a constraint that 
stabilizes pre-critical and post-critical clusters. There is no direct analogy to this method in 
the present theory. The other approach that has more recently begun to be used is steepest- 
descent on the free energy surface. The problem then is that it is necessary to define a 
distance in the parameter space (whether it be the space of density distributions or a space 
of parameters) and there has not been a compelling, much less unique, prescription for doing 
and Q it was proposed that a natural choice was the Euclidean distance 



so. In Refs. 
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in density space, 

^[Pi>P 2 ] = j{ Pl ir)-p 2 {v)) 2 dv (104) 
giving, for arbitrary parameters, 

which should be contrasted with the present results, Eq.( )44|) and Eq.( )68|) . The differences 
between these metrics are not trivial: integration of the steepest-descent equations, Eq.f J76|) . 
using the heuristic g p ab gives completely different results for droplet nucleation from those 
presented above. In fact, the result is the "classical" one whereby a droplet begins as 
a cluster of zero mass and the mass slowly increases as the density and radius increases 
during cluster format ionii^. It is interesting to note that the difference in the results is 
not due to the factor of the inverse density in the expression for the correct metric, but 
rather is due to the fact that it is based on the cumulative mass distribution rather than 
the density. In retrospect, this could be argued to be quite physical: in moving from 
one density configuration to another, what actually is required is the movement of mass 
from one place to another so that the "closeness" of one configuration to another is more 
naturally characterized by how different they are in mass distribution. Furthermore, it is 
to be expected to result from a formulation that respects the conservation of mass and in 
which, therefore, the rate of transport of mass is a limiting factor. Conversely, the incorrect, 
heuristic metric would result from determining the MLP from a dynamics of the form 

^ = -^^ + V^«r,t) (106) 

which is to say, a non-conserved dynamics. (In this equation, Q[p] = F[p] — pN[p] for 
fixed chemical potential, p. In other words, this equation for a non-conserved total mass 
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is to be understood in the grand-canonical ensemble.) While this could be appropriate for 
describing some types of simulations that do not conserve mass, it is clearly inappropriate 
for real physical systems. Finally, it goes without saying that doing steepest-descent using 
an ad-hoc metric based on some parametrization of the free energy is the same thing as 
assuming an ad hoc dynamics which will not, in general, be physically relevant. 

Another point that has not been emphasized is the ability of this formalism to distinguish 
between multiple pathways. Just as there may be multiple local minima, and only one 
true minimum, of a function, there can be multiple local MLP's and only one truly Most 
Likely Path. The simplest and most relevant case occurs when the free energy surface 
possesses multiple saddle points. For example, in the case of protein nucleation there is a 
question of whether the system proceeds directly from low-concentration solution to solid 
or whether it proceeds in two steps by passing through an intermediate, metastable high- 
concentration phase of the type studied above. Until now, the only question that was 
accessible theoretically was the height of the energy barriers and this alone was assumed 
to determine the relevant path. However, this ignores the role of dynamics (for a very 



interesting discussion of this point see Ref. 



53l ) . The present theory allows a way to address 



this problem. First, it is possible to determine local MLP's e.g., in the weak-noise limit 
by means of steepest descent from the various saddle points. Then, given these paths, 
their relative probabilities are easily determined by evaluating Eq. (127j) - (129j) . The absolute 
probability of the paths requires knowledge of possibly unknown normalization constants 
but the relative probabilities can be determined directly by taking the ratio the results. 
This will be elaborated at a later time. 

When applied to globular proteins, namely the case of nucleation of the high-concentration 
amorphous phase from the low-concentration solution, an unexpected result is found whereby 
the process involves two steps. In the first step, a long-wavelength, low-amplitude density 
fluctuation occurs. In the second step, a nucleation event takes place within this region of 
enhanced density. I argued that this is to be expected since, on the one hand, such density 
fluctuations do occur, and at little cost in free energy, and, on the other hand, nucleation 
should be easier in such a region since the free energy barriers should be lower and, equally 
importantly, since the excess mass required for cluster formation is present. This result is 
very robust in that it is independent of the interaction potential and of the free energy- 
functional used while it is valid in the weak-noise limit and could be modified when the 
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effects of strong noise are considered. It is interesting that some of the free-energy-with- 
constraint methods have given hints of a similar pathway, most notably as discussed in Ref. 
and 



12 



16l . It is unclear whether this result is of any practical significance although one might 



speculate that the enhancement or suppression of long-wavelength density (concentration) 
fluctuations could have a strong effect on the nucleation rate. 

The theory developed here is specific to the case of Brownian particles. In some sense, 
however, it can rather be viewed more as an illustration of how a theory can be constructed 
in other cases. For example, for pure fluids, one still has fluctuating hydrodynamics and 
the basic ideas can be developed in a similar way. In that case, temperature fluctuations 
and energy transport are important so that the fluctuating hydrodynamic theory will have 
to be extended to include the heat equation. The main complication might be that the 
hydrodynamic theory cannot be collapsed into a single equation for the density, but rather 
that one must deal will multiple quantities throughout. This may limit the practical utility 
of the theory but it could still yield useful insights into, e.g., the MLP in the weak noise 
limit and the role of the free energy. 

Future directions for further work include the investigation of the role of strong noise. 
Also interesting would be the application to heterogeneous nucleation, nucleation in confined 
systems and the nucleation of solids. 
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Appendix A: Ito and Stratonovich interpretations are the same 

In general, using the Riemann-Stieltjes definition of integration, stochastic integrals differ 
from non-stochastic integrals in that their value depends on where the sampling point is 
chosen within each interval. This means that stochastic differential equations differ according 
to whether they are interpreted in terms of Ito or Stratonovich interpretations. However, 
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in ordinary fluctuating hydrodynamics, the two interpretations turn out to be the same^. 
Here, it is shown that this is also true for the over-damped dynamical model in both the 
general and spherically symmetric cases. Note that this does not hold for the approximate 
order-parameter dynamics thus giving rise to the anomalous force discussed in Section IIII Bl 
and in Appendix [Fl Incidentally, it is worth noting that the fact that the two interpretations 
differ in the latter case is not directly attributable to any approximation: Ito-Stratonovich 
equivalence can be broken simply by a change in variables. 



1. Equivalence for the general model 

In general, the Stratonovich SDE 

-^ = F l (a) + My (a) (t) , (t) (if)) = 5 jt 5 (t - t') (Al) 

corresponds to an Ito SDE of the form 

^ = F> (a) + -^^M^ (a) + M tJ (a) ^ (t) , fe (t) £, (If)) = 8 fl 5 (t - If) (A2) 

so the question comes down to an investigation of the second term on the right in the Ito SDE, 
the so-called spurious drift. To evaluate this for the hydrodynamic model, the procedure 
of van Saarloos^ et al is used whereby the stochastic differential equation is discretized in 
the spatial variable. Specifically, r i — > r n = nA where n = (72.1,712, ■■■,JId) are the (integer) 
coordinates of a lattice point in D dimensions and A is the lattice spacing. Also writing 
p n = p(r n ), the noise term for the over-damped model becomes 

f i 1 D 

V- VTWC M) r -> ^ Yl (VPn+e(")L,n+e(-) ~ y/Pn-&°)£a,n-eM) ( A3 ) 

a=l 

where e^ a ^ is the unit vector in the a-direction (i.e. = <5j a ) and the time-arguments have 
been suppressed. Introducing the notation 5 njI1 / meaning that all components of the vectors 
n and n' are the same the equivalent of My can be identified as 

M n , n / a = — (yPn+eW^n+SKn' ~ VP^SM 5 n-e(«),n') ( A4 ) 



'n+e( a ',n n— e( a ),n 
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We therefore need 



£^f^ (A5) 

kn'a rK 



kn'a * 



k+e( a ),n' ^k-e( a ),n' 



(^n+e( a ),k ~~ ^n-e( a ),k) (^k+e( a ),k ~~ ^k-e( a ),k) 



8A' 

ka 



thus proving equivalence of the two interpretations. 



2. Equivalence for the spherically symmetric model 

In the case of spherical symmetry one has that 

^1 = 4nr*Dp (r) ^^^j + V^W( (r) (A6) 

We discretize as in the previous case but now only require a one-dimensional lattice. Writing 
r i — > r n = nA and keeping in mind that 

4vrr 2 p(r ) = ^M (A7) 
or 

one finds the noise term 



where m n = m(r n ;t), etc. Thus 



^ = ^ -H.-^ (A9) 



o / 7-j vrii-\-\~vrii—\ i 



and so 

(a) = ** ' amf V c """'2a"""' = (A10) 

Appendix B: MLP for potential-driven dynamics with fluctuation-dissipation 
relation 



The purpose of this Appendix is to sketch a straightforward extension of the results of 
Vanden-Eijnden and Heymann^ whereby the assumption of a constant Onsager matrix and 
white noise is lifted while still assuming a fluctuation-dissipation relation. 
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In the following, I consider a set of N stochastic variables, x$ (t) , governed by a diffusive, 
gradient-driven stochastic dynamics with multiplicative noise, 

f = -L (x) ■ (x) + VTecr (x) • £ (f) (Bl) 

where L (x) is a state-dependent matrix of kinetic coefficients, where the scalar constant e 
and matrix a (x) determine the noise amplitude and the noise itself is Gaussian, white and 
diagonally correlated 

<£i(*)£i(0> = M(*-0 (B2) 
Note that the probability density for £ (t) to assume some value, say z, is 

/ X \ N/2 

P (£ (t) = z) = (-) exp (-, 2 /2) (B3) 

The key assumption in the following is that a fluctuation-dissipation relation holds, namely 

L (x) = cr (x) • <r T (x) = D (x) (B4) 

where the second equality reminds that the middle quantity is the diffusion matrix occurring 
in the Fokker-Planck equation. 



Following Refs|56| and |57| the probability density for a given path taking place from time 



t = 0tot = TisP = exp (—^St [x]) , where the action is 

* w = f & + L (x) ' ^ (x) ) ' D " (x) ' (f + L (x) ' (x) ) * (B5) 

and it must be remembered that x depends on time in this and all following expression. 
A path between two points Xj and x 2 is therefore a curve x (t) such that x (0) = Xi and 
x (T) = x 2 . The most likely path (MLP) is determined by minimizing Sq^ [x] over both x (t), 
subject to the constraints on the end points, and time, T. 

In the following, we specialize to the situation that Xi and x 2 are metastable points and 
in fact attractors. We assume that x-space is divided by a separatrix into two regions: region 
I in which points are attracted to xi and region II in which points are attracted to x 2 . 
The separatrix is a curve which will be called S. Any path from xi to x 2 must cross S at 
least once. For the moment, it will be assumed that only one such crossing occurs and the 
possibility of multiple crossings will be discussed below. Any such path can therefore be 
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separated into two pieces: one running from xi to some point x s G S and a second part from 
x s to x 2 : 





[■>'■] 


= S 1 [x 








s 1 


[x] 


-i: 














s 11 













} + S" [x] (B6) 
(f +L (x).^(x)).D-(x).g + L(x).^(x)) df 
5 + L( X ).^ (x) ). D - I(x) .g +L(x) .^ (x) ) di 

with x (T a ) = x s . Clearly, the MLP is determined by minimizing over x s G S and < T s < T 
as well. Consider the second term first. Once the separatrix is crossed all points are attracted 
to x 2 by hypothesis, so the path 

| = -L( X ).^ (X ) (B7) 
x (T s ) = x s 

will eventually reach x 2 : the time required determines T given T s . (Note that one does not 
really start on the separatrix but rather at a point infinitesimally near it on the region II 
side.) However, this path has the property that S n [x] = and this is minimal since the 
integrand is positive definite. This is just the trivial result that the deterministic path is 
the MLP if it passes through the desired points. 

It is not possible to take the deterministic path in region I since the system must go from 
X! to x s and the deterministic dynamics is assumed to always take points in region I towards 
Xi. So, upon noting that expanding the action gives 

S 1 W = T ■ D" 1 (*) ~ + (L (x) • (*)) ■ D- (x) ■ (L (x) ■ i-V (x)) ) dt 

(B8) 



o 



\dt dt \ <9x J \ &k 







Invoking the FDT, the last term is 



K ■ D" 1 (x) ■ L (x) • Ay ( X ) dt = r ^ • ^-V (x) dt (B9) 



dt &k J dt cht 

rTs d 

-^V fx) dt 



dt 

V (x s ) - V (xi) 
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One can therefore write 



S ' H = C 6 ' D " W ' Tt + ( L (x) • S v (x) ) ' D " (x) ' ( L (x) ' l V (x) ) ) * 

(BIO) 

+ 2V (x,)-V(x,) 

= r (f + l (x » ■ <x) ) ■ d_i <x) ■ (i + l <x) ■ i v (x ») * + 4 (x -» - y <xi 

Notice the change in sign of the gradient term. Reversing the sign of the integration variable 
and introducing y (t) — x (— t), so that y (0) = x and y (— T s ) = Xjthe action can be written 

as 

S' N = / (| - L (y) • £v (y)) -D- (y). (I - L (y) . | V (y)) *+4 (V (x.) - V (x,)) 

(Bll) 

The integral is now the action for a path going from xi at t' = —T s to x at t' — so that the 
deterministic path can again be used to set the integral to zero and this will also determine 
T s . One therefore gets that 

S>4(V(x s )- V( Xl )) (B12) 

with equality if the system follows the deterministic paths connecting x s to the end points 
xi and X2. Finally, the action is minimized by choosing x s to be the minimal value on the 
separatrix which is just the critical point. 

Note that recrossing the separatrix will involve a deviation from the deterministic path 
on both sides of the separatrix and so will not minimize the action. For this reason, one 
need only consider a single crossing of the separatrix. 

This serves to establish the claim that for the dynamics given by Eq.( 1Bll) with the 
fluctuation-dissipation relation, Eq. (lB4l) . the MLP crosses the separatrix at the critical point 
and follows the deterministic path away from the critical point. Note the key role played by 
the FDT in this process: the same result will not necessarily hold in the general case. 
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Appendix C: Exact action with order parameters 



Expanding the Lagrangian in Eq. (!52|) gives 



C=- 



2 J r 2 p (r) 



dm (r) 
~~dT 



dr 



(CI) 



—D 
+ 1 2 D '' 



dm (r) d 1 5F [p] 
dt dr r 2 5p (r) 



dr 



r 2 p (r) 



d 1 SF [p] 



dr 



dr r 2 dp (r) 

If the density is parametrized as p(r; t) = p(r; x(t)), then the first term becomes 



2 J r 2 p (r) 



dm (r) 
dt 



dr 



1 



dm (r) dm (r) 



2 Jo r 2 p (r; x) dx a 
. . dx a dxi) 



dx h 



dr 



doc q doc 
~dT~dt 



(C2) 



and the second is 

"°° / 9m(r) d 1 6F[p) 



-D 
-D 



o 



dt dr r 2 5p (r) 
3 d fdm(r) A 5F [p] 

4:71 



dr 



-4tiD^ lim 

dt r-voo 



dt 5p (r) 
9m (r) 5F [p] 



dx a 5p (r) 



dr 



p( r ). 



(C3) 



dr + D 



4irD 



d dm (r] 



dr dt 



4tt 



5F[p] 



5p (r) 



p(r) 



dx a 
~dT 



4irr 



,dp{r) SF[p] 



dx a dp (r) 



dr 



p(r) 



iitD — — +4nD— 

dt dx a dt dx a 

To evaluate this, note that lim r _j.no mfr) = 



iV, the total number of particles in the system. 



Assuming that the boundary condition is that the density assume the bulk value, p^, far 
from the interface gives 



dm (r) 5F [p] 



r'-^-oo \ dx a 5p (r) 



dN 5F[p] 



p( r ). 



dx a 5p(r) 



dN 
dx a 



where the last equality defines the chemical potential. Combining, the final result is 



- D 



dm (r) d 5F [p] 
dt dr 5p (r) 



dr = 4itD 



dx a dQ [p] 
dt dx a 



where Q [p] = F [p] — p,N is the grand potential. Thus, the Lagrangian is 



Air 



-C 



2 9ab[ ^ dt dt 



dx a dxb „ dx a d(3Q [p] 



+ D 



with 



V fx) = -D z 



r 2 p (r) 



dt dx a 



d 1 SF [p] 
dr r 2 5p (r) 



dr 



(C4) 



(C5) 



(C6) 



(C7) 
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Appendix D: Eg. (1801) as a completeness relation 



Assume that {pi (r)}'*L 1 is a complete set of basis functions so that one can write 

oo 

m(r) = ^XiPi (r) (Dl) 



i=l 
-oo 







Xi = m (r) qi (r) dr 
Jo 

where {qi (r)}^ is the bi-orthogonal set satisfying 

Pi ( r ) Qj ( r ) dr = 5ij (D2) 
Then, since m(r)is arbitrary, one has the completeness relation 

S ^2pi{r)qi{r') = 5(r - r') (D3) 

i 

First note that the metric is 

f°° -2 -i , x dm(r)dm(r) , /*oo 

g-y = y r p (r) — ^— fa~ dr = J r P y> P» ( r ) Pi ( r ) rfr ( D4 ) 

The inverse metric is then 

« i =r^ ( °&&*' <d5) 

oo 

r /2 p (r') g.,- (r') g; (r') fir' 
as can be verified by direct evaluation of 

/oo 
r' 2 p(r') qj (r f ) qi (r')dr' (D6) 
j 

/*oo /*oo 

= / dr / dr' (r~ 2 p _1 (r)p,j (r)p.; (r)) (r /2 p (r') gj (r') qi (r')) 

OO /"OO 

dr / dr' (r~ 2 p~ l (r) p^ (r)) 5 (r — r') (r /2 p (r') g^ (r')) 
Jo 



00 



dr pi (r) qi (V) 



ft 



where the second line follows from the completeness relation and the last line by bi- 
orthogonality. Then, one easily verifies that 

? ^^tt* = ^> (r) ( r r " v (r/ ° * (r, ° qj {r " ] dr ") ^ (r ° (d7) 

i 4 3 ii \J / 



"OO 



r //2 p (r") 5 (r - r") 5 (r' - r") dr" 
r 2 p (r) 5 (r — r') 
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as claimed. 

Appendix E: Alternative derivation of order parameter equations 

One way to relate a general functional form, p(r;x(£)), to the exact density, p(r;t), 
is by defining a distance in density space and then minimizing the difference between the 
two functions. It is shown in the main text that the dynamics imposes a unique distance 
functional which is expressed in terms of the corresponding mass functionals as 



Jo 

Determining the fitting parameters at each moment in time by minimizing the difference 
between the actual mass distribution and the approximating distribution gives 



In the following, in order to simplify the notation, the dependence on time will not be 
explicitly indicated. 
Noting that if 




(El) 




(E2) 




(E3) 



then 




(E4) 



so 




(E5) 



where the short-hand notation 




(E6) 



is used. Thus, Eq. flE2p becomes 



= / (m 1 (m (r) ; x) — r) 
Jo 



rrii (m 1 (m (r) ; x) ; x) dm (r) 



dr 



(E7) 



m r (m -1 (m (r) ; x) ; x) dr 
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Another derivative gives 



dt 



00 fhj (m 1 (m (r) 



; x) ; x) rrii (m 1 (m (r) ; x) ; x) dm (r) 



m r (m 1 (m (r) ; x) ; x) m r (m 1 (m (r) ; x) ; x) dr 



+ 



m„ 



m r 



, A/ lit.?- yw» V'^Wl^C 

dm (r) mj (m _1 (m (r) ; x) ; x) dm (r) 



dor,- 



dt 



(fh 1 (m (r) ; x) 



m r \rri~ 1 (m (r) ; x) ; x) dr 

d r^.^-i 
r — 



dt 



(m (m (r) ; x) ; x) 9m (r) 

9r 



m r m 



-1 



m I r ) , xj ; xj 



dr 



(E8) 



or 



' dt 



00 ffl; (fh 1 (m (r) 



,m 



m r 



x) ; x) 9m (r) 



m r (m 1 (m (r) ; x) ; x) m r (m 1 (m (r) 



9r 



dr 



dm (r) m; (m 1 (m (r) ; x) ; x) 9m (r) 



m r 



(m 1 (m (r) ; x) ; x) dt 



dxj 
'"dt" 



(m 1 (m (r) ; x) 



d 

' dt 



:) dr 

rhi (fh~ l (m (r) ; x) ; x) 9m (r) 
(m^ 1 (m (r) , x) ; x) 9r 



dr 



dr 



Expanding in the difference m(r) — m (r; x) gives 



9ij (x (t)) 



dxj (t) 
dt 



9m (r; x (i)) fdm(r;t)\ 1 dm(r;t) 
dxi (t) \ dr J dt 



(E9) 



-dr + O (m (r; t) — m (r; x (t))) 

(E10) 



which gives the result, Eq. (l73|) . 



Appendix F: The spurious drift 



The stochastic model is 



<foi 1 / n 1, , f°° 2D dm (r; x) _ , , , 

= -ZW (x) - e^ 1 (x) / J____lJ-^ (r; t) dr. (Fl) 



dt u 9xj 13 J y47rr 2 p(r;x) 9x.,- 

which must be understood in the Stratonovich interpretation. In order to replace the noise 

term by a simpler form, this must first be written as an Ito SDE. Then, the change of the 

noise can be made and the result transformed back to the Stratonovich interpretation. The 

first step is accomplished using the standard transformation rule^ giving the equivalent Ito 

form of the SDE 

dxi 
~dl 




ID dm (r; x 
Anxr 2 p (r; x) dxj 

dm (r; x 




4irr 2 p (r; x) dxj 
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£ (r; t) dr. 



The autocorrelation of the noise is 

f c 

Da = e 9ij g lm 

Jo 

= ^e 2 gf 



2D dm (r; x) dm (r; x) 



47rr 2 p(r;x) dxj 



dx r 



dr 



(F3) 
(F4) 



so that the noise term in Eq. (lF2p can be replace by the equivalent form ey/ 2Dq i 1 {x) £j(t) 
with g i ;(x)g j i(x) = ^(x) and (f)) = <M(* ~ f)- 

Up to a factor of D, the second term on the right is 



.4 



CO 



dg-/ (x) 

dxj 

+^(x) 9Tk (x) 



(F5) 



1 



47rr 2 p (r; x) \ dxidxi 2 p (r; x) dxi 
Using the symmetry of g^ 1 one has that 

1 d 2 m (r; x) dm (r; x) 



dxi 



9a 1 ( x ) 9ik ( x ) / 
(x) g^ 1 (x) ^ 



47rr 2 p (r; x) dxidx^ dxk 

DO 



(x) fe 1 (x) 



47rr 2 p (r; x) 
1 



d 2 m (r; x) 1 1 <9p (r; x) 9m (r; x) \ dm(r;x) 

dx k 

(F6) 
dr 



-dr. 



-dr 



d 2 m (r; x) dm (r; x) <9 2 m (r; x) <9m (r; x) 



dxidxi 



dx k , 



dxhdxj 



dxi 



d 



2%* ( x ) 9ik (x) 



47rr 2 p (r; x) ctej 
%fc ( x ) 



dxi 



{ dm (r; x) dm (r; x) \ 
V dxi dx k J 
1 <9p (r; x) <9m (r; x) <9m (r; x) 



i.Tir 2 p 2 (r; x) ctej 



dx k 



dr 



so the extra term becomes 



.4 



d9-' 



+ o^/ ( x ) ( x ) 



dxj 

2^( x ) ^ ( x ) 



%fc ( x ) 



(F7) 



/ 4irr 2 f 2 (r; x) 
^ + k'(x)^(x)%M 



<9p (r; x) dm (r; x) <9p (r; x) <9m (r; x) \ dm (r; x 



dxi 



dxi 



dxi 



dxi 



-dr 



dx. 



dxi 



+ 2 (&/ ( x ) 91k 1 ( x ) " ( x ) 9ik ( x )) 



<9p (r; x) dm (r; x) 9m (r; x) 



4 7rr 2p2 jy. x ~j q Xi Q Xk 

Combining these results gives the equivalent Ito form of the SDE with simplified noise term 



dr. 



^ = -2DgJ (x) ^ - De 2 4 T) ~ eV2D g J (x) (t) 



(F8) 



To get the Stratonovich form, one must transform back giving the same equation but with 
spurious drift 



.4 



(S) 



' dx k 



(F9) 
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